******************************************************************************************
******************************************************************************************
************************ BEHAVIORAL RESPONSES  *******************************************

** load data
clear all
use "$maindir/Results/Model_Outputs/ihwap_wedge.dta"

** merge with info about furnace replacements
merge m:1 Household using "$maindir/Data/furn_replace.dta"
drop if _merge==2
drop _merge

***** some extra restrictions on number of homes served by each contractor
gen aux=1
sort Household meterend
bysort Household : gen homeobs = sum(aux) if twoyears_prepost==1 & meterdays_monthly!=. & ProgramYear>=2009
replace homeobs = . if homeobs!=1

sort Household meterend
bysort Household : gen homeobs_post = sum(aux) if treated==1 & twoyears_prepost==1 & meterdays_monthly!=. & ProgramYear>=2009
replace homeobs_post = . if homeobs_post!=1

** fixing builddate to avoid omitted variables
replace builddate = 11 if tab_builddate11==1 | tab_builddate12==1 | tab_builddate13==1

** fixing sqfeet to avoid omitted variables
drop roundsqfeet
gen roundsqfeet = .
	forval i = 0(500)2500 {
		replace roundsqfeet = `i' if sqfeet>`i' & sqfeet<=`i'+500
	}
replace roundsqfeet = 3000 if sqfeet>3000 & sqfeet!=.

** fixing some spending variables to avoid omitted bins
drop binAirCon
gen binAirCon = .
replace binAirCon=0 if Real_tot_actAirCon==0
replace binAirCon=1 if Real_tot_actAirCon>0 & Real_tot_actAirCon!=.

replace binAttic = 2700 if Real_tot_actAttic>2400 & Real_tot_actAttic!=.

** demographic variables
gen simpocc = noccupants
replace simpocc = 5 if noccupants>=5 & noccupants!=.

gen binAge = .
forval i = 20(10)80 {
replace binAge = `i' if Age>`i' & Age<=`i'+10
}

gen binIncome = .
forval i = 5000(5000)35000 {
replace binIncome = `i' if Real_income>`i' & Real_income<=`i'+5000
}
replace binIncome = 1 if Real_income<=5000
replace binIncome = 40000 if Real_income>40000 & Real_income!=.


*** housing structure
gen binSqft = .
forval i = 600(300)2700 {
replace binSqft = `i' if sqfeet>`i' & sqfeet<=`i'+300
}
replace binSqft = 1 if sqfeet<=600
replace binSqft = 3000 if sqfeet>3000 & sqfeet!=.

drop BlowerPreBins
gen BlowerPreBins = .
forval i = 2000(1000)9000 {
replace BlowerPreBins = `i' if Blower_Pre>`i' & Blower_Pre<=`i'+1000
}
replace BlowerPreBins = 1 if Blower_Pre<=2000
replace BlowerPreBins = 10000 if Blower_Pre>10000 & Blower_Pre!=.

*** energy prices
gen binelec = .
local counter = 2
forval i = 10.5(0.5)12.5 {
replace binelec = `counter' if elec_prices>`i' & elec_prices<=`i'+0.5
local counter = `counter' + 1
}
replace binelec = 1 if elec_prices<10.5
replace binelec = `counter' if elec_prices>13 & elec_prices!=.

gen bingas = .
local counter = 2
forval i = 7(1)16 {
replace bingas = `counter' if gas_prices>`i' & gas_prices<=`i'+0.5
local counter = `counter' + 1
}
replace bingas = 1 if gas_prices<7
replace bingas = `counter' if gas_prices>17 & gas_prices!=.

*** transforming decimals in percent
gen savings_gap_pct = savings_gap*100

forval i = 1/200 {
gen savings_gap_pct`i' = 100*savings_gap`i'
}

gen realized_savings_pct = realized_savings*100
winsor2 realized_savings_pct, trim replace cuts(1 99)


forval i = 1/200 {
gen realized_savings_pct`i' = 100*realized_savings`i'
winsor2 realized_savings_pct`i', trim replace cuts(1 99)
}

merge m:1 Household using "$maindir/Data/ihwap_state.dta", keepusing(arch_contID) nogen keep(3)
** keeping number of homes served by contractor
sort arch_contID
by arch_contID : egen totjobs = total(aux) if ww_treat==1
by arch_contID : egen totaljobs = max(totjobs)
drop totjobs
replace totaljobs=. if arch_contID==.
* contractor ID 236 to be ommitted, being one of the biggest contractors in this sample

** main heat pre-treatment
gen MainHeatkBTU = MainHeatBTU/1000

gen binheatkbtu = .
forval i = 60(10)120 {
replace binheatkbtu = `i' if MainHeatkBTU>`i' & MainHeatkBTU<=`i'+10
}
replace binheatkbtu = 1 if MainHeatkBTU<=60
replace binheatkbtu = 130 if MainHeatkBTU>130 & MainHeatkBTU!=.

gen heatworks = 0 if MainHeatOperational=="False"
replace heatworks = 1 if MainHeatOperational=="True"


* fixing variables that have a zero indicator (for later interacting)
foreach x in nstories binAirCon white black hispanic otherrace haselderly hasminor heatworks {
replace `x' = `x'+1
}
replace roundAtticR = 1 if roundAtticR==0

gen heatnotwork = 2 if MainHeatOperational=="False"
replace heatnotwork = 1 if MainHeatOperational=="True"

replace binFurnace = 1 if binFurnace==0

***** some missing vars
gen pct_BlowerReduc = BlowerReduc/Blower_Pre

gen atticR = Attic_RValue
replace atticR = . if atticR==0

gen hasattic = 0 if Attic_RValue==0
replace hasattic = 1 if Attic_RValue>0 & Attic_RValue!=.

winsor2 sqfeet, replace trim cuts(1 99)

gen summer = 0
replace summer = 1 if end_month>=4 & end_month<=9

gen temphdd = HDDbest if treated==0
bysort Household : egen HDDbest_pre = max(temphdd)
drop temphdd
gen temphdd = HDDbest if treated==1
bysort Household : egen HDDbest_post = max(temphdd)
drop temphdd
order HDDbest_pre HDDbest_post, after(HDDbest)

gen HDDbinpre = 1 if HDDbest_pre<50 & HDDbest_pre!=.
replace HDDbinpre = 2 if HDDbest_pre>=50 & HDDbest_pre<55
replace HDDbinpre = 3 if HDDbest_pre>=55 & HDDbest_pre<60
replace HDDbinpre = 4 if HDDbest_pre>=60 & HDDbest_pre<65
replace HDDbinpre = 5 if HDDbest_pre>=65 & HDDbest_pre<70
replace HDDbinpre = 6 if HDDbest_pre>=70 & HDDbest_pre!=.


gen HDDbinpost = 1 if HDDbest_post<50 & HDDbest_post!=.
replace HDDbinpost = 2 if HDDbest_post>=50 & HDDbest_post<55
replace HDDbinpost = 3 if HDDbest_post>=55 & HDDbest_post<60
replace HDDbinpost = 4 if HDDbest_post>=60 & HDDbest_post<65
replace HDDbinpost = 5 if HDDbest_post>=65 & HDDbest_post<70
replace HDDbinpost = 6 if HDDbest_post>=70 & HDDbest_post!=.

gen HDDdiff = HDDbest_pre - HDDbest_post

gen HDDdiffbin = 1 if HDDdiff<-4 & HDDdiff!=.
replace HDDdiffbin = 2 if HDDdiff>=-4 & HDDdiff<-2
replace HDDdiffbin = 3 if HDDdiff>=-2 & HDDdiff<0
replace HDDdiffbin = 4 if HDDdiff==0
replace HDDdiffbin = 5 if HDDdiff>0 & HDDdiff<=2
replace HDDdiffbin = 6 if HDDdiff>2 & HDDdiff<=4
replace HDDdiffbin = 7 if HDDdiff>4 & HDDdiff!=.


*** bins for tmax, tmin, and tavg
gen tavg = (tmax + tmin)/2

gen tmax_f = 1.8*tmax + 32
gen tmin_f = 1.8*tmin + 32
gen tavg_f = 1.8*tavg + 32

gen bintmax = .
forval i = 28(1)89 {
replace bintmax = `i'+1 if tmax_f>`i' & tmax_f<=`i'+1
label define labtmax `=scalar(`i')+1' "(`i' `=scalar(`i')+1']", modify
}
replace bintmax = 28 if tmax_f<=28
replace bintmax = 90 if tmax_f>90 & tmax_f!=.
label define labtmax 28 "<=28", modify
label define labtmax 90 ">90", modify
label values bintmax labtmax

gen bintmin = .
forval i = 10(1)69 {
replace bintmin = `i'+1 if tmin_f>`i' & tmin_f<=`i'+1
label define labtmin `=scalar(`i')+1' "(`i' `=scalar(`i')+1']", modify
}
replace bintmin = 10 if tmin_f<=10
replace bintmin = 70 if tmin_f>70 & tmin_f!=.
label define labtmin 10 "<=10", modify
label define labtmin 70 ">70", modify
label values bintmin labtmin

gen bintavg = .
forval i = 25(1)79 {
replace bintavg = `i'+1 if tavg_f>`i' & tavg_f<=`i'+1
label define labtavg `=scalar(`i')+1' "(`i' `=scalar(`i')+1']", modify
}
replace bintavg = 25 if tavg_f<=25
replace bintavg = 80 if tavg_f>80 & tavg_f!=.
label define labtavg 20 "<=25", modify
label define labtavg 80 ">80", modify
label values bintavg labtavg

* realized savings
gen tau_ml_cv = tau_ml if treated==1
replace tau_ml_cv = cvresids if treated==0
winsor2 tau_ml_cv, replace trim cuts(1 99)

**** alternative binned spending for interactions
gen coarsebinAirSeal = 1 if binAirSeal==1
replace coarsebinAirSeal = 2 if binAirSeal>1 & binAirSeal<=300
replace coarsebinAirSeal = 3 if binAirSeal>300 & binAirSeal<=600
replace coarsebinAirSeal = 4 if binAirSeal>600 & binAirSeal!=.

*** sample to be used in regressions
drop regsample
reghdfe log_total_mmbtu i.treated c.HDD60 c.HDD60#c.HDD60 c.CDD75 c.CDD75#c.CDD75 c.precip c.tmax c.tmin if twoyears_prepost==1 & meterdays_monthly!=. , absorb(i.Household i.month_year) vce(cluster Household end_month)
gen regsample = e(sample) /* to keep stable sample across regressions */


** specification for effect of temperature on NATURAL GAS usage - no controls
reg gas_mmbtu ibn.bintavg if regsample==1 & ProgramYear>=2009 & treated==0 & MainHeatFuel=="1 - Natural Gas", ///
		nocons vce(cluster Household)
est sto tavgbinreg_pre_raw

reg gas_mmbtu ibn.bintavg if regsample==1 & ProgramYear>=2009 & treated==1 & MainHeatFuel=="1 - Natural Gas", ///
		nocons vce(cluster Household)
est sto tavgbinreg_post_raw

coefplot (tavgbinreg_pre_raw, msymbol(Th) msize(*1.2) mcolor(gs4) ciopts(lcolor(gs4))) (tavgbinreg_post_raw, msize(*1.2) mcolor(black) ciopts(lcolor(black))), vertical xtitle("Average Outdoor Temperature ({superscript:o}F)") ytitle("Natural Gas Usage (MMBtu per month)") ///
						graphregion(color(white)) bgcolor(white) ///
						xlabel(5 "(28 29]" 10 "(33 34]" 15 "(38 39]" 20 "(43 44]" 25 "(48 49]" 30 "(53 54]" 35 "(58 59]" ///
						40 "(63 64]" 45 "(68 69]" 50 "(73 74]" 55 "(78 79]", ///
						angle(45) labsize(small)) keep(*bintavg) baselevels ///
						legend(off)	offset(0) ylabel(, glcolor(gray%10))	
graph export "$maindir/Results/Figures/tavgbin_Gas_raw.pdf", replace



********************************************************************************
******************* PRISM ANALYSIS WITH FULL SAMPLE
*** generating multiple HDD bases

local counter = 1
forval i = 55(0.2)65.2 {
	qui : gen hdd_`counter' = `i' - tavg_f
	qui : replace hdd_`counter' = 0 if tavg_f>=`i' & hdd_`counter'!=.
	local counter = `counter'+1
}

matrix hddbases = (.)
forval i = 55(0.2)65.2 {
	sca tempbase = `i'
	matrix hddbases = (hddbases \ tempbase)
}
matrix hddbases = hddbases[2..., 1]


****************** PRISM with controls
matrix prismres_treat_cont = ( . , . , . , . , . )
matrix prismres_notreat_cont = ( . , . , . , . , . )
forval i = 1/51 {

	qui : reg gas_mmbtu hdd_`i' ib900.binFurnace ib300.binAirSeal ib900.binAttic ib1.binWallIns ///
			ib100.binBaseload ib200.binDoor ib1.binWtHtr ib500.binHealSfty ///
			i.binAirCon ib200.binFoundation ib1.binGeneral ib400.binWindow ///
			ib3000.BlowerPreBins ib70.binheatkbtu i.heatworks ///
			ib2.simpocc ib50.binAge ib15000.binIncome ///
			i.white i.black i.haselderly i.hasminor ///
			ib1500.binSqft i.nstories ib1.roundAtticR ib1.builddate ///
			if regsample==1 & ProgramYear>=2009 & treated==0 & MainHeatFuel=="1 - Natural Gas"
	matrix eqn = e(b)
	sca const = eqn[1,2]
	sca slope = eqn[1,1]
	sca rsq = e(r2)
	sca rmse = e(rmse)
	matrix regout = (`=scalar(`i')' , const , slope , rsq , rmse )
	matrix prismres_notreat_cont = (prismres_notreat_cont \ regout)

	qui : reg gas_mmbtu hdd_`i' ib900.binFurnace ib300.binAirSeal ib900.binAttic ib1.binWallIns ///
			ib100.binBaseload ib200.binDoor ib1.binWtHtr ib500.binHealSfty ///
			i.binAirCon ib200.binFoundation ib1.binGeneral ib400.binWindow ///
			ib3000.BlowerPreBins ib70.binheatkbtu i.heatworks ///
			ib2.simpocc ib50.binAge ib15000.binIncome ///
			i.white i.black i.haselderly i.hasminor ///
			ib1500.binSqft i.nstories ib1.roundAtticR ib1.builddate ///
			if regsample==1 & ProgramYear>=2009 & treated==1 & MainHeatFuel=="1 - Natural Gas"
	matrix eqn = e(b)
	sca const = eqn[1,2]
	sca slope = eqn[1,1]
	sca rsq = e(r2)
	sca rmse = e(rmse)
	matrix regout = (`=scalar(`i')' , const , slope , rsq , rmse )
	matrix prismres_treat_cont = (prismres_treat_cont \ regout)
	
	di "HDD base number `i' complete."
}

matrix prismres_notreat_cont = prismres_notreat_cont[2..., 1...]
matrix prismres_treat_cont = prismres_treat_cont[2..., 1...]

matrix rsquares_notreat_cont = prismres_notreat_cont[1..., 4]
matrix rsquares_treat_cont = prismres_treat_cont[1..., 4]

mata : 
rsquares = st_matrix("rsquares_notreat_cont")
i = (.)
w = (.)
maxindex(rsquares, 1, i, w)
st_matrix("maxindex", i)
end
sca maxind_notreat_cont = maxindex[1,1]
sca best_base_notreat_cont = hddbases[`=scalar(maxind_notreat_cont)', 1]

mata : 
rsquares = st_matrix("rsquares_treat_cont")
i = (.)
w = (.)
maxindex(rsquares, 1, i, w)
st_matrix("maxindex", i)
end
sca maxind_treat_cont = maxindex[1,1]
sca best_base_treat_cont = hddbases[`=scalar(maxind_treat_cont)', 1]



***** PRISM with no controls
matrix prismres_treat = ( . , . , . , . , . )
matrix prismres_notreat = ( . , . , . , . , . )
forval i = 1/51 {

	qui: reg gas_mmbtu hdd_`i' if regsample==1 & ProgramYear>=2009 & treated==0 & MainHeatFuel=="1 - Natural Gas"
	matrix eqn = e(b)
	sca const = eqn[1,2]
	sca slope = eqn[1,1]
	sca rsq = e(r2)
	sca rmse = e(rmse)
	matrix regout = (`=scalar(`i')' , const , slope , rsq , rmse )
	matrix prismres_notreat = (prismres_notreat \ regout)

	qui: reg gas_mmbtu hdd_`i' if regsample==1 & ProgramYear>=2009 & treated==1 & MainHeatFuel=="1 - Natural Gas"
	matrix eqn = e(b)
	sca const = eqn[1,2]
	sca slope = eqn[1,1]
	sca rsq = e(r2)
	sca rmse = e(rmse)
	matrix regout = (`=scalar(`i')' , const , slope , rsq , rmse )
	matrix prismres_treat = (prismres_treat \ regout)
	
	di "HDD base number `i' complete."
}

matrix prismres_notreat = prismres_notreat[2..., 1...]
matrix prismres_treat = prismres_treat[2..., 1...]

matrix rsquares_notreat = prismres_notreat[1..., 4]
matrix rsquares_treat = prismres_treat[1..., 4]

mata : 
rsquares = st_matrix("rsquares_notreat")
i = (.)
w = (.)
maxindex(rsquares, 1, i, w)
st_matrix("maxindex", i)
end
sca maxind_notreat = maxindex[1,1]
sca best_base_notreat = hddbases[`=scalar(maxind_notreat)', 1]

mata : 
rsquares = st_matrix("rsquares_treat")
i = (.)
w = (.)
maxindex(rsquares, 1, i, w)
st_matrix("maxindex", i)
end
sca maxind_treat = maxindex[1,1]
sca best_base_treat = hddbases[`=scalar(maxind_treat)', 1]



*****************************************************************************
******** plotting PRISM rsquared versus HDD bases
tempfile ihwapprism
save `ihwapprism'

**** with controls
clear
svmat prismres_notreat_cont
svmat prismres_treat_cont

rename (prismres_notreat_cont1 prismres_notreat_cont2 prismres_notreat_cont3 prismres_notreat_cont4 prismres_notreat_cont5) ///
		(Base Constant Beta R_sq RMSE)

rename (prismres_treat_cont1 prismres_treat_cont2 prismres_treat_cont3 prismres_treat_cont4 prismres_treat_cont5) ///
		(Base2 Constant2 Beta2 R_sq2 RMSE2)

local counter = 1
forval i = 55(0.2)65.2 {
local hddbase : di %9.1f `i'
label define baselab `counter' "`hddbase'", modify
local counter = `counter' + 1
}
label values Base baselab
label values Base2 baselab

sum R_sq, detail
sca maxrsq = r(max)
sum Base if R_sq==maxrsq
sca hddbest = r(mean)
sum R_sq2, detail
sca maxrsq2 = r(max)
sum Base2 if R_sq2==maxrsq2
sca hddbest2 = r(mean)

twoway (scatter R_sq Base, msymbol(Th) msize(*1.2) mcolor(gs8))  ///
	(scatter R_sq2 Base2, msize(*1.2) mcolor(black)) ///
	(pci `=scalar(maxrsq)' 0 `=scalar(maxrsq)' `=scalar(hddbest)', lcolor(gs8)) ///
	(pci 0.71 `=scalar(hddbest)'  `=scalar(maxrsq)' `=scalar(hddbest)', lcolor(gs8)) ///
	(pci `=scalar(maxrsq2)' 0 `=scalar(maxrsq2)' `=scalar(hddbest2)', lcolor(black)) ///
	(pci 0.71 `=scalar(hddbest2)'  `=scalar(maxrsq2)' `=scalar(hddbest2)', lcolor(black)), ///
	xlabel(1(20)51, valuelabel labsize(small) angle(90)) ///
	xlabel(`=scalar(hddbest)', add valuelabel) ///
	ylabel(`=round(maxrsq, 0.0001)', add angle(0)) ///
	xlabel(`=scalar(hddbest2)', add valuelabel) ///
	ylabel(`=round(maxrsq2, 0.0001)', add angle(0)) ///
	xtitle("Heating Degree Day Base") ytitle("PRISM R-Squared") ///
	legend(off) graphregion(color(white)) bgcolor(white) ylabel(, glcolor(gray%10))
graph export "$maindir/Results/Figures/hddbest_prepost_controls.pdf", replace 


**** no controls
clear
svmat prismres_notreat 
svmat prismres_treat

rename (prismres_notreat1 prismres_notreat2 prismres_notreat3 prismres_notreat4 prismres_notreat5) ///
		(Base Constant Beta R_sq RMSE)

rename (prismres_treat1 prismres_treat2 prismres_treat3 prismres_treat4 prismres_treat5) ///
		(Base2 Constant2 Beta2 R_sq2 RMSE2)

local counter = 1
forval i = 55(0.2)65.2 {
local hddbase : di %9.1f `i'
label define baselab `counter' "`hddbase'", modify
local counter = `counter' + 1
}
label values Base baselab
label values Base2 baselab

sum R_sq, detail
sca maxrsq = r(max)
sum Base if R_sq==maxrsq
sca hddbest = r(mean)
sum R_sq2, detail
sca maxrsq2 = r(max)
sum Base2 if R_sq2==maxrsq2
sca hddbest2 = r(mean)

twoway (scatter R_sq Base, msymbol(Th) msize(*1.2) mcolor(gs8))  ///
	(scatter R_sq2 Base2, msize(*1.2) mcolor(black)) ///
	(pci `=scalar(maxrsq)' 0 `=scalar(maxrsq)' `=scalar(hddbest)', lcolor(gs8)) ///
	(pci 0.64 `=scalar(hddbest)'  `=scalar(maxrsq)' `=scalar(hddbest)', lcolor(gs8)) ///
	(pci `=scalar(maxrsq2)' 0 `=scalar(maxrsq2)' `=scalar(hddbest2)', lcolor(black)) ///
	(pci 0.64 `=scalar(hddbest2)'  `=scalar(maxrsq2)' `=scalar(hddbest2)', lcolor(black)), ///
	xlabel(1(20)51, valuelabel labsize(small) angle(90)) ///
	xlabel(`=scalar(hddbest)', add valuelabel) ///
	ylabel(`=round(maxrsq, 0.0001)', add angle(0) labsize(small)) ///
	xlabel(`=scalar(hddbest2)', add valuelabel) ///
	ylabel(`=round(maxrsq2, 0.0001)', add angle(0)) ///
	xtitle("Heating Degree Day Base") ytitle("PRISM R-Squared") ///
	graphregion(color(white)) bgcolor(white) ylabel(, glcolor(gray%10)) ///
	legend(order(1 "Pre-Treatment" 2 "Post-Treatment"))
graph export "$maindir/Results/Figures/hddbest_prepost.png", replace width(2500)
graph export "$maindir/Results/Figures/hddbest_prepost.pdf", replace 


clear
use `ihwapprism'

** model with optimal HDD base - pre-treatment
qui: gen hdd_simul = hdd_`=scalar(maxind_notreat_cont)'
qui: reg gas_mmbtu hdd_simul ib900.binFurnace ib300.binAirSeal ib900.binAttic ib1.binWallIns ///
		ib100.binBaseload ib200.binDoor ib1.binWtHtr ib500.binHealSfty ///
		i.binAirCon ib200.binFoundation ib1.binGeneral ib400.binWindow ///
		ib3000.BlowerPreBins ib70.binheatkbtu i.heatworks ///
		ib2.simpocc ib50.binAge ib15000.binIncome ///
		i.white i.black i.haselderly i.hasminor ///
		ib1500.binSqft i.nstories ib1.roundAtticR ib1.builddate ///
		if regsample==1 & ProgramYear>=2009 & treated==0 & MainHeatFuel=="1 - Natural Gas"
drop hdd_simul
est sto tempmod_pre
matrix tempcoefs_pre = e(b)
	
** model with optimal HDD base - post-treatment
qui: gen hdd_simul = hdd_`=scalar(maxind_treat_cont)'
qui: reg gas_mmbtu hdd_simul ib900.binFurnace ib300.binAirSeal ib900.binAttic ib1.binWallIns ///
		ib100.binBaseload ib200.binDoor ib1.binWtHtr ib500.binHealSfty ///
		i.binAirCon ib200.binFoundation ib1.binGeneral ib400.binWindow ///
		ib3000.BlowerPreBins ib70.binheatkbtu i.heatworks ///
		ib2.simpocc ib50.binAge ib15000.binIncome ///
		i.white i.black i.haselderly i.hasminor ///
		ib1500.binSqft i.nstories ib1.roundAtticR ib1.builddate ///
		if regsample==1 & ProgramYear>=2009 & treated==1 & MainHeatFuel=="1 - Natural Gas"
drop hdd_simul
est sto tempmod_post
matrix tempcoefs_post = e(b)



********************************************************************************
*************** simulate impact of changing only the slope (efficiency)

** counterfactual predictions using pre-treatment PRISM model
gen hdd_simul = hdd_`=scalar(maxind_notreat_cont)'
est restore tempmod_pre
predict temppred_counterfactual if regsample==1 & ProgramYear>=2009 & treated==1 & MainHeatFuel=="1 - Natural Gas"
drop hdd_simul

** post-treatment predictions according to PRISM model
gen hdd_simul = hdd_`=scalar(maxind_treat_cont)'
est restore tempmod_post
predict temppred_post if regsample==1 & ProgramYear>=2009 & treated==1 & MainHeatFuel=="1 - Natural Gas"
drop hdd_simul

** predictions using slope from post, but HDD and baseload from pre-treatment PRISM model
gen hdd_simul = hdd_`=scalar(maxind_notreat_cont)'
est restore tempmod_pre
matrix b = tempcoefs_pre
matrix b[1,1] = tempcoefs_post[1,1]
ereturn post b
predict temppred_newslope if regsample==1 & ProgramYear>=2009 & treated==1 & MainHeatFuel=="1 - Natural Gas"
drop hdd_simul

** savings from the slope change
sum temppred_newslope if regsample==1 & ProgramYear>=2009 & treated==1 & MainHeatFuel=="1 - Natural Gas" & savings_gap!=.
sca avgpred_newslope = r(mean)

sum temppred_counterfactual if regsample==1 & ProgramYear>=2009 & treated==1 & MainHeatFuel=="1 - Natural Gas" & savings_gap!=.
sca avgpred_counterfactual = r(mean)

sca avgsave_newslope = avgpred_newslope - avgpred_counterfactual
sca avgpctsave_newslope = avgsave_newslope/avgpred_counterfactual

** total savings accoring to PRISM approach
sum temppred_post if regsample==1 & ProgramYear>=2009 & treated==1 & MainHeatFuel=="1 - Natural Gas" & savings_gap!=.
sca avgpred_post = r(mean)

sca avgsave_post = avgpred_post - avgpred_counterfactual
sca avgpctsave_post = avgsave_post/avgpred_counterfactual

** share of savings attributed to slope change
sca slopeshare = avgsave_newslope/avgsave_post

di "Average MMBtu savings due to slope change = `=round(scalar(avgsave_newslope), 0.01)'"
di "Average percent savings due to slope change = `=round(scalar(100*avgpctsave_newslope), 0.01)'%"
di "Average MMBtu savings according to PRISM approach = `=round(scalar(avgsave_post), 0.01)'"
di "Average percent savings according to PRISM approach = `=round(scalar(100*avgpctsave_post), 0.01)'%"
di "Share of savings attributed to the slope change = `=round(scalar(100*slopeshare), 0.01)'%"

matrix appendixf2 = (`=round(scalar(avgsave_newslope), 0.01)' , `=round(scalar(100*avgpctsave_newslope), 0.01)' , `=round(scalar(avgsave_post), 0.01)' , `=round(scalar(100*avgpctsave_post), 0.01)' , `=round(scalar(100*slopeshare), 0.01)')
matrix rownames appendixf2 = "Value"

esttab matrix(appendixf2) using "$maindir/Results/Tables/appendixf2.tex", ///
	mlabels(, none) collabels("Average MMBtu savings due to slope change" "Average percent savings due to slope change" "Average MMBtu savings according to PRISM approach" "Average percent savings according to PRISM approach" "Share of savings attributed to the slope change") ///
	replace 


********************************************************************************
*************** simulate potential impact of behavior on energy usage

** no controls
gen hdd_simul = hdd_`=scalar(maxind_treat)'
reg gas_mmbtu hdd_simul if regsample==1 & ProgramYear>=2009 & treated==1 & MainHeatFuel=="1 - Natural Gas"
predict pred_base_`=scalar(maxind_treat)'
drop hdd_simul

forval i = 1(1)6 {
gen hdd_simul = hdd_`=scalar(maxind_treat) - `i''
predict pred_base_`=scalar(maxind_treat) - `i''
drop hdd_simul
}


** with controls
gen hdd_simul = hdd_`=scalar(maxind_treat_cont)'
reg gas_mmbtu hdd_simul ib900.binFurnace ib300.binAirSeal ib900.binAttic ib1.binWallIns ///
		ib100.binBaseload ib200.binDoor ib1.binWtHtr ib500.binHealSfty ///
		i.binAirCon ib200.binFoundation ib1.binGeneral ib400.binWindow ///
		ib3000.BlowerPreBins ib70.binheatkbtu i.heatworks ///
		ib2.simpocc ib50.binAge ib15000.binIncome ///
		i.white i.black i.haselderly i.hasminor ///
		ib1500.binSqft i.nstories ib1.roundAtticR ib1.builddate ///
		if regsample==1 & ProgramYear>=2009 & treated==1 & MainHeatFuel=="1 - Natural Gas"
gen controlsamp = e(sample)		
predict pred_cont_base_`=scalar(maxind_treat_cont)'
drop hdd_simul

forval i = 1(1)6 {
gen hdd_simul = hdd_`=scalar(maxind_treat_cont) - `i''
predict pred_cont_base_`=scalar(maxind_treat_cont) - `i''
drop hdd_simul
}


****** summarizing the results in energy and gap space
* with controls
forval i = 0(1)6 {
gen total_mmbtu_cont_base`=scalar(maxind_treat_cont) - `i'' = electric_mmbtu + pred_cont_base_`=scalar(maxind_treat_cont) - `i''
gen tau_ml_base`=scalar(maxind_treat_cont) - `i'' = total_mmbtu_cont_base`=scalar(maxind_treat_cont) - `i'' - pred_full
gen realized_savings_base`=scalar(maxind_treat_cont) - `i'' = tau_ml_base`=scalar(maxind_treat_cont) - `i''/pred_full if treated==1
gen savings_gap_cont_base`=scalar(maxind_treat_cont) - `i'' = realized_savings_base`=scalar(maxind_treat_cont) - `i'' - projected_savings
winsor2 savings_gap_cont_base`=scalar(maxind_treat_cont) - `i'', replace trim cuts(0.5 99.5)
}

foreach x in  savings_gap savings_gap_cont_base35 savings_gap_cont_base34 savings_gap_cont_base33 savings_gap_cont_base32 savings_gap_cont_base31 savings_gap_cont_base30 savings_gap_cont_base29 {
	reg `x' i.aux if regsample==1 & ProgramYear>=2009 & treated==1 & MainHeatFuel=="1 - Natural Gas" & savings_gap_cont_base35!=. & savings_gap!=., nocons 
est sto est`x'
}

foreach x in  realized_savings realized_savings_base35 realized_savings_base34 realized_savings_base33 realized_savings_base32 realized_savings_base31 realized_savings_base30 realized_savings_base29 {
	reg `x' i.aux if regsample==1 & ProgramYear>=2009 & treated==1 & MainHeatFuel=="1 - Natural Gas" & savings_gap_cont_base35!=. & savings_gap!=., nocons 
est sto est`x'
}

drop pred_base_35-savings_gap_cont_base29 


******************************************************************
******* bootstrapping for PRISM analysis

*** Mata function that is used within loop - this is necessary since calling Mata breaks a loop in Stata
mata : 
void findmax(real scalar treat) {
	real matrix i
	real matrix w
	real matrix rsquares
	
	if (treat == 0) {
		rsquares = st_matrix("rsquares_notreat_cont")
	}
	else {
		rsquares = st_matrix("rsquares_treat_cont")
	}
	
	i = (.)
	w = (.)
	maxindex(rsquares, 1, i, w)
	st_matrix("maxindex", i)
	}
end


matrix fullbootsmat_gap = (., ., ., ., ., ., ., .)
matrix fullbootsmat_save = (., ., ., ., ., ., ., .)
forval j = 1/200 {

	di "starting bootstrap iteration `j'"
		
	****************** PRISM with controls
	matrix prismres_treat_cont = ( . , . , . , . , . )
	matrix prismres_notreat_cont = ( . , . , . , . , . )
	forval i = 21/50 {
		di "            starting HDD base `i'"

		qui : reg gas_mmbtu hdd_`i' ib900.binFurnace ib300.binAirSeal ib900.binAttic ib1.binWallIns ///
				ib100.binBaseload ib200.binDoor ib1.binWtHtr ib500.binHealSfty ///
				i.binAirCon ib200.binFoundation ib1.binGeneral ib400.binWindow ///
				ib3000.BlowerPreBins ib70.binheatkbtu i.heatworks ///
				ib2.simpocc ib50.binAge ib15000.binIncome ///
				i.white i.black i.haselderly i.hasminor ///
				ib1500.binSqft i.nstories ib1.roundAtticR ib1.builddate ///
				[fweight = weight_boots`j'] if regsample==1 & ProgramYear>=2009 & treated==0 & MainHeatFuel=="1 - Natural Gas"
		matrix eqn = e(b)
		sca const = eqn[1,2]
		sca slope = eqn[1,1]
		sca rsq = e(r2)
		sca rmse = e(rmse)
		matrix regout = (`=scalar(`i')' , const , slope , rsq , rmse )
		matrix prismres_notreat_cont = (prismres_notreat_cont \ regout)

		qui : reg gas_mmbtu hdd_`i' ib900.binFurnace ib300.binAirSeal ib900.binAttic ib1.binWallIns ///
				ib100.binBaseload ib200.binDoor ib1.binWtHtr ib500.binHealSfty ///
				i.binAirCon ib200.binFoundation ib1.binGeneral ib400.binWindow ///
				ib3000.BlowerPreBins ib70.binheatkbtu i.heatworks ///
				ib2.simpocc ib50.binAge ib15000.binIncome ///
				i.white i.black i.haselderly i.hasminor ///
				ib1500.binSqft i.nstories ib1.roundAtticR ib1.builddate ///
				[fweight = weight_boots`j'] if regsample==1 & ProgramYear>=2009 & treated==1 & MainHeatFuel=="1 - Natural Gas"
		matrix eqn = e(b)
		sca const = eqn[1,2]
		sca slope = eqn[1,1]
		sca rsq = e(r2)
		sca rmse = e(rmse)
		matrix regout = (`=scalar(`i')' , const , slope , rsq , rmse )
		matrix prismres_treat_cont = (prismres_treat_cont \ regout)

	}

	matrix prismres_notreat_cont = prismres_notreat_cont[2..., 1...]
	matrix prismres_treat_cont = prismres_treat_cont[2..., 1...]

	matrix rsquares_notreat_cont = prismres_notreat_cont[1..., 4]
	matrix rsquares_treat_cont = prismres_treat_cont[1..., 4]

	mata : findmax(0)
	
	sca maxind_notreat_cont = maxindex[1,1]
	sca best_base_notreat_cont = hddbases[`=scalar(maxind_notreat_cont)', 1]

	mata : findmax(1)
	
	sca maxind_treat_cont = maxindex[1,1]
	sca best_base_treat_cont = hddbases[`=scalar(maxind_treat_cont)', 1]




	*********** simulations
	** with controls
	qui: gen hdd_simul = hdd_`=scalar(maxind_treat_cont)'
	qui: reg gas_mmbtu hdd_simul ib900.binFurnace ib300.binAirSeal ib900.binAttic ib1.binWallIns ///
			ib100.binBaseload ib200.binDoor ib1.binWtHtr ib500.binHealSfty ///
			i.binAirCon ib200.binFoundation ib1.binGeneral ib400.binWindow ///
			ib3000.BlowerPreBins ib70.binheatkbtu i.heatworks ///
			ib2.simpocc ib50.binAge ib15000.binIncome ///
			i.white i.black i.haselderly i.hasminor ///
			ib1500.binSqft i.nstories ib1.roundAtticR ib1.builddate ///
			[fweight = weight_boots`j'] if regsample==1 & ProgramYear>=2009 & treated==1 & MainHeatFuel=="1 - Natural Gas"
	qui: predict pred_cont_base_`=scalar(maxind_treat_cont)'
	drop hdd_simul

	forval i = 1(1)6 {
	qui: gen hdd_simul = hdd_`=scalar(maxind_treat_cont) - `i''
	qui: predict pred_cont_base_`=scalar(maxind_treat_cont) - `i''
	drop hdd_simul
	}


	****** summarizing the results in gap space
	* with controls
	forval i = 0(1)6 {
		gen total_mmbtu_cont_base`=scalar(maxind_treat_cont) - `i'' = electric_mmbtu + pred_cont_base_`=scalar(maxind_treat_cont) - `i''
		gen tau_ml_base`=scalar(maxind_treat_cont) - `i'' = total_mmbtu_cont_base`=scalar(maxind_treat_cont) - `i'' - pred_boots`j'
		gen realized_savings_base`=scalar(maxind_treat_cont) - `i'' = tau_ml_base`=scalar(maxind_treat_cont) - `i''/pred_boots`j' if treated==1
		gen savings_gap_cont_base`=scalar(maxind_treat_cont) - `i'' = realized_savings_base`=scalar(maxind_treat_cont) - `i'' - projected_savings
		winsor2 savings_gap_cont_base`=scalar(maxind_treat_cont) - `i'', replace trim cuts(0.5 99.5)
	}

	matrix coefmat_gap = (.)
	foreach x in savings_gap savings_gap_cont_base`=scalar(maxind_treat_cont)' savings_gap_cont_base`=scalar(maxind_treat_cont) - 1' savings_gap_cont_base`=scalar(maxind_treat_cont) - 2' ///
			savings_gap_cont_base`=scalar(maxind_treat_cont) - 3' savings_gap_cont_base`=scalar(maxind_treat_cont) - 4' ///
			savings_gap_cont_base`=scalar(maxind_treat_cont) - 5' savings_gap_cont_base`=scalar(maxind_treat_cont) - 6' {
		qui: reg `x' i.aux [fweight = weight_boots`j'] if regsample==1 & ProgramYear>=2009 & treated==1 & MainHeatFuel=="1 - Natural Gas" & savings_gap_cont_base`=scalar(maxind_treat_cont)'!=. & savings_gap!=., nocons 
		matrix tempmat = e(b)
		matrix coefmat_gap = (coefmat_gap, tempmat)
	}
	matrix coefmat_gap = coefmat_gap[1..., 2...]
	matrix fullbootsmat_gap = (fullbootsmat_gap \ coefmat_gap)
	
	matrix coefmat_save = (.)
	foreach x in realized_savings realized_savings_base`=scalar(maxind_treat_cont)' realized_savings_base`=scalar(maxind_treat_cont) - 1' realized_savings_base`=scalar(maxind_treat_cont) - 2' ///
			realized_savings_base`=scalar(maxind_treat_cont) - 3' realized_savings_base`=scalar(maxind_treat_cont) - 4' ///
			realized_savings_base`=scalar(maxind_treat_cont) - 5' realized_savings_base`=scalar(maxind_treat_cont) - 6' {
		qui: reg `x' i.aux [fweight = weight_boots`j'] if regsample==1 & ProgramYear>=2009 & treated==1 & MainHeatFuel=="1 - Natural Gas" & savings_gap_cont_base`=scalar(maxind_treat_cont)'!=. & savings_gap!=., nocons 
		matrix tempmat = e(b)
		matrix coefmat_save = (coefmat_save, tempmat)
	}
	matrix coefmat_save = coefmat_save[1..., 2...]
	matrix fullbootsmat_save = (fullbootsmat_save \ coefmat_save)
	
	drop pred_cont_base_`=scalar(maxind_treat_cont)'-savings_gap_cont_base`=scalar(maxind_treat_cont) - 6' 

}
matrix fullbootsmat_gap = fullbootsmat_gap[2..., 1...]
matrix fullbootsmat_gap = 100*fullbootsmat_gap

matrix fullbootsmat_save = fullbootsmat_save[2..., 1...]
matrix fullbootsmat_save = 100*fullbootsmat_save


matrix meanmat = fullbootsmat_gap[1..., 2]
matrix ppmat = fullbootsmat_gap[1..., 2...]
mata : 
A = st_matrix("meanmat")
B = st_matrix("ppmat")
st_matrix("ppmat", mm_colvar(B))
C = B :- A
D = C :/ A
D = 100*D
D
st_matrix("percentmat", mm_colvar(D))
end

do "$maindir/Code/bootstrapres.do"
forval i = 0(1)6 {
	local colnum = `i' + 1
	local basenum = 35 - `i'
	matrix tempmat = ppmat[1, `colnum']
	matrix bootsvar`basenum'  = tempmat
	mat colnames bootsvar`basenum' = "c1"
	mat rownames bootsvar`basenum' = "c1"
	
	matrix tempmat = percentmat[1, `colnum']
	matrix bootspctvar`basenum'  = tempmat
	mat colnames bootspctvar`basenum' = "c1"
	mat rownames bootspctvar`basenum' = "c1"
	
	est restore estsavings_gap_cont_base`basenum'
	matrix tempmat = e(b)
	scalar tempscalar = tempmat[1,1]
	matrix bootscoef`basenum' = (100*tempscalar)
	sca nobs = e(N)
	mat colnames bootscoef`basenum' = "c1"
	mat rownames bootscoef`basenum' = "c1"
	scalar scalar`basenum' = bootscoef`basenum'[1,1]
	
	matrix bootspctcoef`basenum' = 100*(scalar`basenum' - scalar35)/scalar35
	mat colnames bootspctcoef`basenum' = "c1"
	mat rownames bootspctcoef`basenum' = "c1"
	
	bootstrapres bootscoef`basenum' bootsvar`basenum' nobs
	est sto resbase`basenum'
	
	bootstrapres bootspctcoef`basenum' bootspctvar`basenum' nobs
	est sto pctresbase`basenum'
}


matrix meanmat = fullbootsmat_save[1..., 2]
matrix ppmat = fullbootsmat_save[1..., 2...]
mata : 
A = st_matrix("meanmat")
B = st_matrix("ppmat")
st_matrix("ppmat", mm_colvar(B))
C = B :- A
D = C :/ A
D = 100*D
D
st_matrix("percentmat", mm_colvar(D))
end


forval i = 0(1)6 {
	local colnum = `i' + 1
	local basenum = 35 - `i'
	matrix tempmat = ppmat[1, `colnum']
	matrix bootsvar`basenum'  = tempmat
	mat colnames bootsvar`basenum' = "c1"
	mat rownames bootsvar`basenum' = "c1"
	
	matrix tempmat = percentmat[1, `colnum']
	matrix bootspctvar`basenum'  = tempmat
	mat colnames bootspctvar`basenum' = "c1"
	mat rownames bootspctvar`basenum' = "c1"
	
	est restore estrealized_savings_base`basenum'
	matrix tempmat = e(b)
	scalar tempscalar = tempmat[1,1]
	matrix bootscoef`basenum' = (100*tempscalar)
	sca nobs = e(N)
	mat colnames bootscoef`basenum' = "c1"
	mat rownames bootscoef`basenum' = "c1"
	scalar scalar`basenum' = bootscoef`basenum'[1,1]
	
	matrix bootspctcoef`basenum' = 100*(scalar`basenum' - scalar35)/scalar35
	mat colnames bootspctcoef`basenum' = "c1"
	mat rownames bootspctcoef`basenum' = "c1"
	
	bootstrapres bootscoef`basenum' bootsvar`basenum' nobs
	est sto savebase`basenum'
	
	bootstrapres bootspctcoef`basenum' bootspctvar`basenum' nobs
	est sto pctsavebase`basenum'
}


esttab resbase35 resbase34 resbase33 resbase32 resbase31 resbase30 resbase29 ///
				using "$maindir/Results/Tables/new_rebound_results.tex", ///
				replace b(3) se(3) nostar nonumbers ///
				mlabels("61.8" "61.6" "61.4" "61.2" "61.0" "60.8" "60.6") ///
				nonotes stats(N, fmt(%9.0gc) labels("Observations"))

esttab pctresbase35 pctresbase34 pctresbase33 pctresbase32 pctresbase31 pctresbase30 pctresbase29 ///
				using "$maindir/Results/Tables/new_rebound_pctresults.tex", ///
				replace b(3) se(3) nostar nonumbers ///
				mlabels("61.8" "61.6" "61.4" "61.2" "61.0" "60.8" "60.6") ///
				nonotes stats(N, fmt(%9.0gc) labels("Observations"))

esttab savebase35 savebase34 savebase33 savebase32 savebase31 savebase30 savebase29 ///
				using "$maindir/Results/Tables/new_reboundsave_results.tex", ///
				replace b(3) se(3) nostar nonumbers ///
				mlabels("61.8" "61.6" "61.4" "61.2" "61.0" "60.8" "60.6") ///
				nonotes stats(N, fmt(%9.0gc) labels("Observations"))

esttab pctsavebase35 pctsavebase34 pctsavebase33 pctsavebase32 pctsavebase31 pctsavebase30 pctsavebase29 ///
				using "$maindir/Results/Tables/new_reboundsave_pctresults.tex", ///
				replace b(3) se(3) nostar nonumbers ///
				mlabels("61.8" "61.6" "61.4" "61.2" "61.0" "60.8" "60.6") ///
				nonotes stats(N, fmt(%9.0gc) labels("Observations"))

